%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% Explanation of what the file does
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% This function solves for sectoral prices (hats or dots) and level trade 
% shares in the RUV economy with forward looking migration using a 
% contraction mapping algorithm
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% Inputs needed
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% 1. Initial guess for prices (startguess), size I*S
% 2. Previous trade shares (PastSha), size I*I*S
% 3. Dot or hat trade costs (TAUmatrix), size I*I*S
% 4. Dot or hat technologies (TECmatrix), size I*S
% 5. Dot or hat wages (wagessec),  size I*S
% 6. Labor shares (labshares), size I*S
% 7. IO matrix (inoutmat), size I*S*S
% 8. Trade elasticities (Smatrix), size 1*S
% 9. Max number of iterations (maxiter), scalar
% 10. Tolerance for stopping (tolerance), scalar
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% Outputs produced
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% 1. Dot or hat prices (NPrice), size I*S
% 2. New trade shares in levels (NTradeShare), size I*I*S

function [NPrice,NTradeShare]=NEBA_PriceSolver(startguess,PastSha,...
TAUmatrix,TECmatrix,wagessec,labshares,inoutmat,Smatrix,maxiter,tolerance)

% First, define some important variables

pnewguess      = startguess;
I              = size(pnewguess,1);
S              = size(pnewguess,2);
iterations     = 1;
distancebg     = 10;

% Second, run the contraction mapping loop for sectoral prices (hat or dot)

while distancebg  > tolerance && iterations<maxiter
    pguess     = pnewguess;
    interpart  = squeeze(prod(repmat(pguess,1,1,S).^inoutmat,2));
    priceown   = TECmatrix .* (wagessec.^labshares .* interpart)...
                 .^(1-repmat(Smatrix,I,1));
    costtemp   = TAUmatrix.^(1-repmat(permute(Smatrix,[3,1,2]),I,I,1))...
                 .* repmat(permute(priceown,[1,3,2]),1,I,1);
    pnewguess  = (squeeze(sum(PastSha.*costtemp,1)))...
                 .^(1./(1-repmat(Smatrix,I,1)));
    iterations = iterations+1;
    distancebg = max(abs(pnewguess(:)-pguess(:)));
end

% Third, assign final price matrix to the output and define trade shares

NPrice         = pnewguess;
NTradeShare    = (PastSha.*costtemp)./repmat(sum(PastSha.*costtemp),I,1,1);

% Fourth, print a warning if iterations are close to or at limit

if iterations > maxiter-3
    disp("Maximum number of iterations in price contraction reached")
end

end
